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Abstract 

We present results of Monte-Carlo simulations for finite 2D single and bilayer systems. 
Strong Coulomb correlations lead to arrangement of particles in configurations resembling a 
crystal lattice. For binary layers, there exists a particularly rich variety of lattice symmetries 
which depend on the interlayer separation d. We demonstrate that in these mesoscopic 
lattices there exist two fundamental types of ordering: radial and orientational. The 
dependence of the melting temperature on d is analyzed, and a stabilization of the crystal 
compared to a single layer is found. 

Introduction. The properties of a finite number of charged particles (N ~ 10) in a 
single "two-dimensional" layer have been the subject of intensive threoretical and ex- 
perimental investigation in the last decade 0-0. In particular, if at low temperature 
the density parameter r s is increased, transitions from a Fermi liquid to a so-called 
Wigner molecule [|J and, further, to a Wigner crystal have been found ||. The inter- 
play between the long-range Coulomb interaction and a shallow confinement potential 
plays in these systems a prominent role making conventional effective single-particle 
approximations unreliable. In this connection the predictions of both classical and 
quantum Monte Carlo methods which treate N-body correlations rigorously, are of 
great importance, in particular for the theoretical understanding of the solid-liquid 
crossover. Recently, a new 2D system has attracted the attention of several groups, 
namely, bilayer structures. The phase diagram of bilayer systems is far more rich 
compared to single layers. In particular, the formation of Wigner lattices in electronic 
or ionic bilayers has been predicted, both by classical and quantum-mechanical @ 
studies and revealed the existence of distinct structural phases: rectangular, square, 
rhombic and triangular staggered lattices. The transition between these phases takes 
place at specific values of the interlayer distance d, when one of the phases becomes 
energetically favorable and may proceed in continuous or discontinuous manner. How- 
ever, these predictions have been made for macroscopic systems. Here, we extend the 
analysis of bilayer stuctures to finite-size systems (e.g. electrons or ions in quantum 
dots and radio frequency traps, respectively). Our previous inverstigation of finite sin- 
gle layer structures || has revealed that these mesoscopic systems have a richer phase 
diagram than their macroscopic counterpart, and that the phase diagram strongly 
depends on the particle number. In this paper, we demonstrate that different crys- 
talline structures are stable in different ranges of the interlayer distance. Moreover, 
we find that melting may proceed in several stages, and the crystal may be stabilized 
by choosing a proper value of d. 

Model and characteristic parameters. We consider a system of N charged particles 
of the same type interacting via the repulsive Coulomb potential located in two 2D 
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layers which are a distance d apart. In each layer a circular symmetric parabolic 
potential of strength uj is applied. The system is described by the hamiltonian 

N t2 v 2 N m *,, 2 r 2 2 Nl p 2 1 Nl N2 p 2 

i=l 4 i=l 1=1 i<j 01 lJI i=l j=l \/ 3 ^ 

where r^- = rj — r^, m* and e& are the effective electron mass and background dielectric 
constant, respectively. We use the following length and energy scales: r , given by 
e 2 / e 6 r o = m*u 2 r 2 t /2, and E c - the average Coulomb energy, E c = e 2 /eb r o- After the 
scaling transformations {r — > r/ro, — > E/E c , d — *> d/ro} the hamiltonian becomes 

2 N AT 2 V, JVi JV 2 

« = yE v i+E r . 2 + EEuq + ^EE^=- w 

i=i <=i ;=i i<,- 1 i=i j=i Jrjj + d 2 

where n = v / 2^o/ r o = ( a B/ r o) ; a s is the effective Bohr radius, and Zjj = h/(m*uj ), 
is the extension of the ground state wave function of noninteracting trapped electrons. 
Further, we define, in analogy to macroscopic systems, r s = ro/a* B = 1/n 2 0. Finally, 
we introduce the dimensionless temperature T = k B T / E c which allows us to define 
the classical coupling parameter as V = l/T 0. To make reliable calculations in the 
crystal phase, where the Coulomb energy strongly exceeds the kinetic energy, stan- 
dard methods, such as Hartree-Fock or density functional theory are not applicable. 
We, therefore, use classical and path integral (PIMC) Monte Carlo simulations in the 
classical and quantum regions of the phase diagram, respectively. 

Phase boundary of the mesoscopic Wigner crystal in a single layer, d = 0. Our 
simulations revealed that the structure of the clusters in the crystal phase strongly 
depends on the particle number and results from the competition of two ordering ten- 
dencies: for large particle numbers, a triangular lattice is energetically favorable. In 
contrast, for small N, the particles tend to form shells. For intermediate particle num- 
bers, iV > 40, only the outer electrons form shells, while the inner part goes over to a 
triangular lattice structure. These configurations are remarkably stable and are visible 
even outside the crystal phase. To distinguish the solid and liquid phases, we used 
the standard Lindemann criterion which is based on the analysis of the magnitude of 
the inter-particle distance fluctuations, see e.g. 0. Indeed, our simulations show that 
these fluctuations exhibit jumps (see Fig|3].a below) which allow us to locate the phase 
boundary of the crystal. When leaving the crystal, the relative distance fluctuations 
increase several times which results from particles exchanging their lattice sites or un- 
dergoing inter-shell transitions. Interestingly, we find two distinguished crystal phases: 
first, a completely ordered state and second, a partially (only radially) ordered phase 
where shells can rotate with respect to each other. The transition between the two 
phases will be called orientational melting ("OM"), whereas the transition from the 
radially ordered state to the liquid-like state is radial melting ("RM"). We underline 
that this phase still resembles the crystal (it still exhibits shell structure). Neverthe- 
less, it is characterized by essentially increased inter-shell exchanges of particles and, 
therefore, this "Wigner molecule" phase 0] may not be identified with a crystal. 
The results for the phase boundaries in the density-temperature plane are summarized 
in Fig. |l] for various particle numbers. Consider first the line of radial melting "RM" . 
At low densities, n < 0.03, we are in the classical regime, and the phase boundary 
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is given by critical values of the coupling parameter T which, in Fig. [I] corresponds 
to straight lines from the origin. We found that F deviates from the 2D bulk value, 
Too = 137, e.g. JJ], and strongly depends on the cluster size. Further increase of 
T, (reduction of the slope of the straight line), eventually leads to the second phase 
boundary, labelled "OM" which corresponds to freezing of the inter-shell rotation. At 
sufficiently low temperature, the electron behavior in this phase is dominated by quan- 
tum effects. This becomes particularly clear if the density parameter n is increased. 
In this case, we observe increasing quantum (zero-point) fluctuations. When the line 
"OM" is reached, again a jump of the relative distance fluctuations is observed which 
corresponds to "cold" orientational melting ||. Finally, increasing n further, leads to 
the line "RM", i.e. to cold radial melting of the crystal. From the figure it is clear 
that each of the two crystal phases exists up to a maximum temperature. This tem- 
perature is comparatively low. Using typical parameters of 2D GaAs semiconductors, 
the maximum temperature of the radially ordered phase is in the range of 1 . . . 5K. 




Fig. 1: Phase diagram of the mesoscopic 2D Wigner 
crystal. "OM" ("RM") denotes the orientational 
(radial) melting curves for N = 12, 19, 20. Insert 
shows an enlarged picture of the low-density region. 
Dotted straight lines indicate the radial melting 
transition of a macroscopic classical and quantum 
WC. Brueckner parameter follows from the density 
by r s = 1/n 2 . Shown error bars are typical for all 



Interestingly, the phase bound- 
ary of the crystal significantly varies 
with the number of particles. This 
dependence is particularly strong for 
the orientational melting transition, 
and we find that crystals in so-called 
magic clusters have unusually high 
stability. The reason is that these 
clusters have the highest angular 
symmetry (the particle numbers on 
the shells have a common divisor). 
This is a peculiarity of mesoscopie 
systems; the strong number depen- 
dence of the melting properties van- 
ishes in the limit of a macroscopic 
system where also the two crystal 
phases merge into one. 



curves. 



Bilayers, d ^ 0. Let us now 
turn to electron crystallization in 
two parallel 2D layers. We study 
a classical bilayer cystal using stan- 
dard MC simulations. The particles 
are arranged into two parallel layers each containing an equal number of particles. For 
zero separation, d — 0, we have just one 2D "atom". For example, for 2N = 38, this 
is a cluster with one electron in the trap center and three shells containing 7, 13 and 
17 electrons, respectively. In the opposite limit, d —>■ oo, the system consists of two 
independent clusters with N = 19 particles and shell configuration {1,6,12}. Below 
we provide a detailed investigation of the stuctural changes between these two limits 
for the case N = 19. 



Different crystalline structures can be identified from snapshots of instantaneous 
configurations (see insets in Fig. |2|). In Fig. |3|.c we also show particle configurations 
in both layers projected onto the same plane. To characterize the stuctural symmetry 
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in each layer we use a suitable order parameter given by 



G ra (i?) = ^X)|0m(i)ly 
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Fig. 2: Order parameter G m (R), Eq. (||), characterizing 
the angular lattice symmetry for several values of d. 
Symbols denote m = 4 (diamonds), m = 5 (squares) 
and m = 6 (triangles). Insets show the corresponding 
particle configuration in one layer. 



where the first sum is taken over 
the Ni particles in the inner part 
of the crystal (all particles ex- 
cept those on the outer shell), 
and the sum over k runs over 
all intralayer neighbors of par- 
ticle j within a circle of ra- 
dius R. m = 2,3...; and 
6jk is the angle between some 
fixed axis and the vector con- 
necting the j — th particle and 
its k—t\i neighbor. For a trian- 
gular lattice, the order parame- 
ter G e approaches one, whereas 
for a square (rhombic) lattice, G 4 
becomes close to unity. We 
also consider the magnitude of 
the relative intra/inter- layer dis- 
tance fluctuations. In the 
vicinity of orientational and ra- 
dial melting, the fluctuations, 
Eq. (|]), show a strong increase, 
thus providing a suitable quan- 



titative criterion for these phase transitions (see Fig. [|a,b). We define the interlayer 
(/ 7^ m, C = NiN m ) and intralayer [I = m and i ^ j, C = Ni(Ni — 1)] distance 



fluctuations, 



u 



Ini 



and the radial fluctuations with respect to the trap center, u , as 



where iVj and N m are the number of particles in layer / and m, respectively; is the 
projection of the distance between particles i and j onto one of the layers; and (. . .) 
denotes an ensemble average. 

First, we start to decrease the interlayer distance from the limit d — *■ oo. The 
configuration {1,6,12} corresponding to two independent clusters does not change 
up to (I k 0.9. The temperature dependence of the inter(intra)-layer pair distance 
fluctuations u lm (u l ), Eq. (f|), for d = 1.4 are shown in Fig. [|a. The two jumps 
in the behavior of the fluctuations correspond to the two melting transitions. First, 
orientational melting takes place in both layers simultaneously. That means that the 
shells with 6 and 12 particles are orientationally disordered and can rotate relative to 
each other, as in the single-layer system. This melting takes place at a temperature 
T fa 2.5 ■ 10~ 3 . Notice that the intra and interlayer fluctuations differ significantly 
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which shows that interlayer correlations are comparatively small which allows the shells 
in the two layers to rotate independently of each other (Fig. [3|.a, upper part). Radial 
melting sets in at a higher temperature, T r pa 6.5 • 10~ 3 , which is the same temperature 
as was found for single-layer systems p. 

Decreasing the in- 
terlayer distance to d — 
0.9 leads to redistri- 
bution of particles on 
the shells. The cluster 
configuration {1,7,11} 
in each layer becomes 
energetically favorable, 
cf. Fig. |2|. These 
changes in the clus- 
ter symmetry lead to 
quite different temper- 
atures of orientational 
melting. Now the inner 
shells in the two layers 
are "frozen" and can- 
not rotate with respect 
to each other. How- 
ever, they can rotate to- 
gether relative to the 
(frozen) two outer shells 
which is clearly seen 
in the coinciding inter 
and intralayer distance 
fluctuations, Fig. |[a, 
lower part. This is a 
new type of disordering 
transition which is miss- 
ing in the single-layer 

system. The critical temperature of this transition is significantly lower, T pa 5 • 10~ 6 , 
which is a result of the modified shell configuration. Total melting occures at 
T r pa 7 ■ 10" 3 . 

Now we consider the distance d = 0.8. The interlayer correlations lead to the for- 
mation of staggered rhombic lattices in the inner parts of the clusters. The interesting 
point here is that these staggered lattices still have the possibility to rotate relative 
to the outer shells. This takes place at T Q pa 4.5 ■ 10~ 3 . This essential increase of the 
orientational melting temperature (compared to the case d = 0.9) is due to the fact 
that the outer and inner shells begin to lose angular symmetry (which is similar to 
d = 0.5, see inset of Fig. ^]). Close inspection of this figure shows that some particles 
(in the top part) are moved inbetween the shells which effectively hampers rotation of 
the shells. The calculated radial distribution functions (not shown) also confirm this 
conclusion. In addition to the peaks corresponding to the two shells (for d > 0.9), here 
a third maximum arises inbetween. Radial melting sets in at T r pa 7.5 • 10 -3 , which 
is an important increase of the stability against radial disordering compared to larger 
values of d, cf. Fig. |3].b. 



Fig. 3: Crystallization in a classical bilayer system for varying 
interlayer distance d. Left three figures: Relative two-particle 
distance fluctuations u lm for particles from the same layer (open 
symbols) and from different layers (full symbols) and radial fluc- 
tuations u\ Eq. (4). Upper right figure shows snapshots of the 
crystal structure for T = 3 • 10~ 3 . Full and open symbols denote 
particles from different layers, thin lines are guide for the eye 
to underline the cluster symmetry. Lower right figure shows the 
critical temperature of the radial (RM) and orientational (OM) 
melting transitions versus d. 
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At interparticle distances in the range d = 0.5 . . . 0.7 the angular shell symmetry 
in the system is practically lost and the particle configuration in each layer resembles 
a rhombic lattice, cf. Fig. |3].c. This leads to drastic changes in the behavior of the 
fluctuations: orientational and radial transitions take place approximatly at the same 
temperature, T Q « 7 . . . 8 • 10~ 3 and T r 7.5 . . . 9 • 10~ 3 , cf. Fig. |||.d. Further decrease 
of the interlayer separation up to d = 0.3 gives rise to another symmetry change, to 
formation of a rectangular lattice in each layer, cf. Fig. |2|. Interestingly however, the 
system as a whole regains spherical symmetry: as is clearly seen in Fig. [3|.c, all particles 
together arrange into a spherically symmetric cluster of 2N = 38 particles with the 
shell configuration {1,7,13,17}. This restoration of spherical symmetry causes the 
temperature of the orientational melting to drop significantly around d = 0.3, cf. 
Fig||d. The obvious conclusion is that the interlayer correlations already dominate the 
behavior of the clusters, and the system has become effectively a single-layer structure. 
When d — > 0, we again observe all properties of a single-layer crystal, including the 
critical values of temperatures for radial and orientational melting, T Q ~ 1 • 10~ 6 and 
T r « 7.2 • 10~ 3 , respectively. 

In summary, we have presented a detailed analysis of Wigner crystallization of 
finite electron systems in one and two layers. In particular, we discussed the influence 
of interlayer correlations on the crystal phase. We have shown that, when d is reduced, 
the clusters first lose spherical symmetry and transform into a rectangular (rhombic) 
lattice. Finally, at low d, a spherical arrangement is re-established when all particles 
effectively form a single layer. We have found that this intermediate range of <i-values 
provides an additional stabilization of the Wigner crystal, and the melting temperature 
may rise by up to 50% (confirming the previous finding for macroscopic systems ||). 
The same tendency is expected to hold in mesoscopic clusters in the quantum range 
(very low temperature): interlayer correlations should significantly reduce the critical 
value of r s || up to r s — > 2/3r s which should improve prospects for an experimental 
observation of electron (hole) crystallization in semiconductor heterostructures even at 
zero magnetic field. 
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